Pseudorandom number generators and the square site percolation threshold 
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A select collection of pseudorandom number generators is applied to a Monte Carlo study of the 
two dimensional square site percolation model. A generator suitable for high precision calculations is 
identified from an application specific test of randomness. After extended computation and analysis, 
an ostensibly reliable value of Pc = 0.59274598(4) is obtained for the percolation threshold. 
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I. INTRODUCTION 

The square site percolation threshold, Pc, is a clearly 
and simply defined mathematical concept [l|, 0]- Per- 
colation models have been well studied, and are known 
for their numerous applications pi]. Yet to date, no 
analytical expression has been found for the numerical 
value of Pc- The square site lattice lacks the symme- 
try that has allowed exact solutions on other topologies 
[2, [j, M, M, 0, H ly] ■ So long as the problem remains in- 
tractable, statistical estimates from Monte Carlo studies 
can, at least, offer approximate values. Such calculations 
invariably make extensive use of some form of pseudo- 
random number generator (PRNG). 

A PRNG is a deterministic algorithm that outputs 
a sequence of words with properties closely mimick- 
ing those of a truly random sequence. Well analysed 
generators include the linear congruential, lagged Fi- 
bonacci, generalised feedback shift register, and deriva- 
tives thereof Jp, [H [H, [H, [11 M, M O M, M M 
[2ll . [22 . [23l . [24| . Because these algorithms are simple, 
they do not produce output with the complexity of a ran- 
dom sequence [23, l20|. The autocorrelation coefficients 
of a pseudorandom sequence are not identically zero, and 
these departures from true randomness introduce a sam- 
pling bias that leads to systematic error. 

Twenty years ago, concern was being given to the de- 
mands then made of PRNGs in calculations using 10^^ 
pseudorandom numbers generated at MHz rates [23] ■ Re- 
cently, high performance parallel computer systems with 
thousands, rather than tens or hundreds, of processors 
have become much more widely available. These enable 
calculations with 10^^ pseudorandom numbers generated 
at GHz rates, and are likely to play a central role in fu- 
ture research. Very high precision can now be achieved 
through brute force of sampling, but accuracy is another 
matter. For reliable Monte Carlo estimates at these new 
higher precision levels, the PRNG(s) chosen must be of 
sufficient quality. Hence contemporary demands upon 
PRNGs are, and will continue to become, much greater 
than in the past. 

This study compares several established PRNGs 
within the context of the square site percolation problem. 
Following application specific testing, a seemingly reli- 
able generator is identified. This is subsequently used to 
locate the percolation threshold with, in principle, both 



accuracy and precision. 



II. GENERATORS 

Throughout this study, the computational word 
length, w, shall be fixed at 32. All arithmetical oper- 
ations taking place within any PRNG are performed in 
modulo 2"". All PRNG arithmetical operands, and prod- 
ucts thereof, are members of {0:2"" — 1}, where {a : b} 
denotes the set of all integers not less than a and not 
greater than b. Consequently the words of any PRNG 
output sequence also belong to {0 : 2™ — 1}. The «th 
word of an output sequence shall be denoted by Xi . With 
one noted exception, no output sequence is decimated in 
any way. The first million words of each sequence are 
discarded prior to beginning any Monte Carlo sampling 
procedure. 

Some PRNGs make use of bit-wise operations within 
their internal algorithms. The notation adopted here is 
© for bit-wise Boolean logical exclusive-or, and >m for 
shift m bits to the right (where m is a positive integer). 
Whenever these bit-wise operations are performed, the 
operands are decomposed into their respective standard 
binary representations, most-significant bit (leftmost) to 
least-significant bit (rightmost). Arithmetic being con- 
strained to a subset of the integers, any bits shifted to a 
position right of the decimal point are lost. Hence, within 
this study, the operation of Ottj is equivalent to integer 
division by 2™. 

The specific PRNGs considered within this exercise are 
defined as follows. 

TT is the two-tap additive lagged Fibonacci genera- 
tor Xi = Xi-4is + Xi-\279- This generator has previously 
been used for high-precision percolation threshold mea- 
surement by Newman and Ziff [2^, 123] ■ 

TTT combines the output from a pair of two-tap gen- 
eralised feedback shift-register generators, Ui = Ui-^-ji ® 
Wi-9689 and Vi — Wi-30 © 'Ui-127, to return a single word 
Xi = Ui © Vi. This is the generator most likely used for 
two and three dimensional percolation by Deng and Blote 

usim. 

SWB is a Marsaglia and Zaman subtract with borrow 
generator, Xi = Xi-222 — Xi-237 — /3i-i, where the bor- 
row. Pi, is equal to one if Xi-222 < 2:1-237 + A-i, and is 
otherwise equal to zero [la, 123] • 



QTA is the quad-tap generalised feedback shift-register 
generator Xi = a;i_|57®a;i_3i4©a;i_47i ©Si-ggsg, as used 
by Ziff and Stell [11 [^ (see [13). 

QTB is the quad-tap generahsed feedback shift-register 
generator Xi = a;i_47i ©Xi_i586 ® Xi_6988 ® a^i-gesQ- This 
generator has been used by Newman and Ziff, and has 
been found to produce threshold estimates consistent 
with those of the TT generator p^. l29l. l3i| . 

XG is Brent's xorgen4096 generator [23] ■ Specifically, 
the implementation xorgen4096i, from his C language 
xorgens304 distribution, was that used here. This gener- 
ator has performed well in randomness tests conducted 
by L'Ecuyer and Simard [24| . 

MT is Matsumoto and Nishimura's MT19937 
Mersenne twister generator [1^. Specifically, their 
MT19937ar C language distribution was the implemen- 
tation used here. The MT19937 algorithm has been 
used for computing integrals in semi-ri gor ous work by 
Balister, Bollobas, Walters and Riordan [3^, [S^], and for 
Monte Carlo sampling by Lee [37l |. 

DMT is a pair of MT generators operated entirely in- 
dependently of one another. The output sequence from 
each of these generators is decimated, with only every 
fourth word used. Lattice sites are then selected by 
means of their Cartesian coordinates, using one number 
from each generator. This scheme has previously been 
used by Lee [33]. 

Let L be the number of sites lying on each edge of 
an L X i square lattice. Every site on that lattice is 
typically given a unique index or label, j 6 {0 : L^ — 1}. 
In a microcanonical ensemble Monte Carlo calculation 
[23, [2^, |3J] ' such as those performed here, PRNG output 
words, Xi, are used to pseudorandomly select sites, Sj, for 
occupation. Now, consider a transformation T{x, N) = 
X [> {w — log2 N) , where A'^ is a positive-integer power 
of two. This is a distribution preserving many-to-one 
surjective map from the integers x £ {0 : 2™ — 1} to 
the integers T{x,N) G {0 : N — 1}. Further consider a 
bijection H that maps the integers {0 : A^ — 1} onto site 
labels. With the usual choice of site labels also being 
the integers {0 : N — 1}, H is conventionally taken to 
be the identity map. In the single generator systems 
defined above, PRNG output word Xi is associated with 
site Sj(^^.) via j{xi) = H{T{xi,L'^)). For the DMT, a 
pair of PRNG output words, Ui and Vi (one from each 
of the output decimated MT generators), is mapped to 
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H{T{u„L)+LT{v,,L)). This 



halves the number of bits actually used from each output 
word (the most significant bits being those retained). 

Each of these generators must be provided with an 
initial finite sequence of words from which to begin cal- 
culating an infinite pseudorandom sequence. In the case 
of the TT generator for instance, a list of some 1279 
initial words is required. These initial lists were con- 
structed by one of four simpler generators, here denoted 
LCGa, LCGb, LCGm and WMx. LCGa is the hnear 
congruentialgenerator, Xi = 69069xi_i -I- 1, suggested by 
Marsaglia [12i] . LCGb is a similar linear congruential gen- 



erator, Xi = 69069a::i_i -I- 1234567, also due to Marsaglia 
j2(j ]. LCGm is the modified linear congruential genera- 
tor, Xi = 1812433253(xi_i © {xi^i [> 30)) + i, appearing 
in Matsumoto and Nishimura's MT19937ar distribution 
of their Mersenne twister algorithm. WMx is the Weyl 
modified Marsaglia xorshift generator built into the xor- 
gens4096i algorithm appearing within Brent's xorgens304 
distribution [21, 23]. These initialisation generators are 
themselves seeded from a single word, xq G {0 : 2™ — 1}. 
The TTT and DMT generators both require two initial- 
isation lists, each being derived from one of these four 
generators, each starting with a distinct independent seed 
word. 



III. TEST PROCEDURE 

The above listed generators were compared, in the con- 
text of site percolation on the square lattice, by using 
each one to make a Monte Carlo estimate of the crossing 
probability function, Rl^u, a,t L — 2048, over the do- 
main n G {2474000 : 2498300}. i?L,„ is defined as the 
probability that a single cluster connects two specified 
opposing boundary sides of the N — L x L square lattice 
in the microcanonical ensemble when precisely n random 
sites are occupied. The value of i?2048.n monotonically 
increases from around 0.05 at n = 2474000 to around 
0.95 at n = 2498300. Hence the occupation domain 
studied encompasses the critical region of the percolative 
phase transition. The numerous lattice configurations 
required to accurately determine RL,n were constructed, 
from each PRNG output sequence, over the above do- 
main only, by the unbiased algorithm of Lee [371 ]. The 
only exception was the XG generator from which samples 
were obtained over the same domain by the unbiased al- 
gorithm of Newman and ZiflF [2^, [2^ . 

The Newman and Ziff binomial convolution 



N 



RLip) = J2( P"(l-p)'^""i?L. 



(1) 



then gives the crossing probability, Rl{p), in the canoni- 
cal ensemble where each lattice site is independently ran- 
domly occupied at probability p |28l l29l |34| . In principle 
the summation should run over all n G {0 : N}, but 
as samples were taken only over the restricted domain 
of n above, the summation was truncated accordingly. 
The standard deviation of the binomial distribution in 
equation ([T]) is given by OL.p ~ L^p{l — p). The data 
analysis here is concerned with values of p such that the 
distribution maximum, located at n — nint(pA^), lies be- 
tween 10 and 12 (yL,p from the nearest end of the sam- 
pling region. Consequently, the truncation induced error 
in Rl{p) is not more than 10^^^. This is completely 
negligible when compared to statistical sampling uncer- 
tainties, which were never less than 10~®. 

The canonical crossing probability curve is used to 
identify a site occupation probability, pf (L), defined such 



that 



i?i(pf(L)) = l/2 + fc/L, 



(2) 



where k — 0.320(1), as determined by Ziff and Newman 
[3J| (their parameter bo). For p « Pc to first order 
Rl{p) ~ 0.5 + k/L + {{p - Pc)L^/'') [33, and hence 
P({L) provides a reasonable estimate of the critical point 
Pc. Ziff and Newman have found that the second order 
equation i?L(pc) ~ 0.5+fcL~^ — 0.44_L^^ is abetter model 
of the data at small L [3^, however the L~^ term is neg- 
ligible for L > 1024 at the levels of precision considered 
here. Values for pf (2048) were thus obtained from each 
of the PRNGs described above. These were subsequently 
compared against each other and against previous Pc es- 
timates made with the same generators. For large L, 
Rl{p) rises very steeply in the neighbourhood oi p k, p^. 
Consequently, pi{L) is relatively insensitive to the exact 
value of k provided that k/L <^ 1. The uncertainty in k 
limits the maximum attainable precision in pf(2048) to 
±2 X 10^9. 

Combinatorial terms of the binomial distribution in 
equation ([1]) were calculated by the essentially exact 
method of Newman and Ziff [29|. For p near pc, use 
of the Gaussian approximation to the binomial would 
have introduced an error of order 10~® in i?2048(p)j this 
corresponding to an error of order 10^ ^'^ in p itself. It 
is sometimes possible to dispense with the convolution 
altogether and make a microcanonical ensemble approxi- 
mation oiRLip = n/N) w RL,n- With L = 2048, and for 
p near p^, this introduces an error of around 4 x 10^^ in 
Rl{p): which corresponds to an error of around 2 x 10^^ 
in p. This approximation is acceptable at low enough 
precision, has the advantage that only a much narrower 
domain of sampling need be considered, and has been 
employed in earlier work by Lee [37|. However, since 
the induced error (measured as the difference between 
p and n/N such that either Rl{p) — Rl.u — 0.5 or 
Rl{p) = Rl.u = 0.5 -I- k/L) was found to scale as only 
^-1.5(3)^ when a set of measurements arc to be taken over 
a range of lattice sizes, to precisions of order 1/N, the 
error will become significant at large L. The approxima- 
tion was not adopted here. 

Correlations inevitably found in the output sequence of 
any deterministic pseudorandom number generator will 
result in correlations within the spatial pattern of occu- 
pied sites upon the lattice. As noted by Compagner [3q . 
this in turn will bias the resulting Monte Carlo estimate 
of the crossing probability function. Consequently, when 
estimates obtained from two different generators are in- 
consistent, then at least one of those generators likely suf- 
fers from significant correlations in its output sequence, 
hence rendering it unsuitable for use at the level of preci- 
sion of the study. Because the true values of Rl.u, Rl{p) 
and Pc are unknown, it will be unclear as to which of the 
generators is deficient. 

While there is merit in performing general tests on 
PRNGs, it is often preferable to have an application spe- 
cific test such as the sensitive hull walk of Ziff [13] . Here 



FIG. 1: Standard (left) and (example) nonstandard (right) 
enumerations of the L — 4 square lattice. 

So si S2 S3 So sii se si 

54 S5 Sq S7 Si2 S7 S2 S13 

55 Sg Sio Sii Ss S3 Si4 Sg 
S12 Si3 Si4 Si5 Si Si5 Sio S5 



a scheme is used that changes the relation between nu- 
merical PRNG output sequences and the spatial patterns 
of occupied lattice sites, without altering the underlying 
problem or topology in any way. The standard enumer- 
ation of the lattice, as shown in figure [1] prescribes a 
specific relation between patterns in PRNG output and 
in clusters of occupied sites. By adopting some other 
(nonstandard) enumeration, as per the example in figure 
[1] some different relation is obtained. An ideal random 
number generator will produce results independent of the 
chosen enumeration. A pseudorandom number genera- 
tor, with correlations in its output sequence, will produce 
results that do depend upon the enumeration. By com- 
paring results from a common generator on two different 
lattice enumerations, inadequate, outcome biasing, gen- 
erators may be identified. This simple application spe- 
cific test does not require knowledge of the percolation 
threshold or spanning probability curves. 

The direct approach to implementing such an enumer- 
ation is to allocate each site a set of pointers explicitly 
identifying its geometrical neighbours. In the nonstan- 
dard enumeration of figure (TJ for example, S7 would have 
pointers to S2, S3, su and S12. Sites are pseudorandomly 
selected as per normal and the Monte Carlo sampling 
proceeds just as for the standard enumeration. In prac- 
tice this results in a dramatic performance decrease of the 
simulations (more than a factor of two was found in this 
study) . The problem is believed to be the cache prefetch- 
ing of the high performance computer system used, where 
if, in some linear array, Sj is being accessed then the hard- 
ware assumes that Sj+i (being the next contiguous data 
element in memory) will be wanted next. 

An alternative method is to change the mapping, H, 
between scaled PRNG output words, y — T{x,N), and 
site labels, j. In the standard enumeration of figure [I] 
H(jj) = y is the identity map. In the nonstandard enu- 
meration, H{0) = 0, H{11) = 1, H{6) ^ 2, and so on, 
with the general relation being H{y) — 3y (mod 16). 
This is analogous to the cluster label labels of Hoshen 
and Kopelman (39[. The nonstandard enumeration is 
that effectively in use, while the standard enumeration 
is preserved in computer memory, thus avoiding perfor- 
mance problems. 

When the hash function, H, is simple (that is, of sim- 
ilar algebraic complexity to the PRNG), it will not so 
much hide correlations in the output sequence as man- 
ifest those patterns in some other way, giving rise to a 
different estimate for p{{L). If, on the (systematic) stan- 
dard enumeration, correlations in PRNG output lead to 
spatial correlations of occupied sites that in turn bias the 



estimate, then, on some other (systematic) nonstandard 
enumeration, those same PRNG correlations will give rise 
to spatial correlations of a different nature that bias the 
estimate in some other way. This provides a simple, ap- 
plication specific, test for PRNG biasing of the Monte 
Carlo samples. If a given PRNG is correlation free, then 
the estimates derived from it will be independent of the 
lattice enumeration. If instead, the PRNG output does 
suffer from correlations, then different enumerations may 
lead to different results. The test can be tuned, with the 
hash chosen so as to maximise the observed shift in the 
test quantity. A simple hash related to the taps or period 
is bound to highlight intrinsic PRNG shortcomings [3q . 
Alternatively, a hash much more complex than the gener- 
ator algorithm could go some way toward hiding output 
sequence correlation induced bias. 

Hence define two further generators, TTH and MTH, 
as (respectively) the exact same TT and MT generators 
defined previously, but with a somewhat arbitrary non- 
trivial mapping H{y) = 947j/ (mod N) between integers 
y e {0 : iV - 1} and site labels j e {0 : N - 1}. On a 
lattice oi L — 2048, this hash is equivalent to a system- 
atic nonstandard enumeration where the rightward and 
downward neighbours of site Sj are (when they exist) 
Sj-i-3171195 (mod N) and Sj+i824768 (mod N) respectively. 



IV. TEST RESULTS 

The various estimates of pf(2048) thus obtained are 
listed in table HI Results are separated according to the 
generator and initialisation scheme used. SWBb, for in- 
stance, indicates the SWB generator with its initialisa- 
tion list derived from the LCGb output sequence. Simi- 
larly, XGx is the XG generator initialised from the WMx 
output sequence. Results shown are based on surveys of 
order 10^ effectively independent samples per occupation 
level n. Each of these sets involved the generation of or- 
der 10^^ to 10^^ (approaching 10^^ in the xorgens case), 
pseudorandom numbers. 

TTa is the TT generator initialised from LGGa. The 
TTa based pf (2048) estimate in table U is consistent with 
the results of Newman and Ziff that were also obtained 
(primarily [40|) from the TT generator (see table HI)) . 
TTHa is the hashed TT generator, again initialised with 
LGGa. The TTa and TTHa based estimates are suffi- 
ciently different to indicate the probable existence of sta- 
tistically significant correlations within the TT generator 
output sequence. This gives cause for concern about the 
use of the TT PRNG for this application at this level of 
precision. 

TTTab is the TTT generator with its initial u and 
V lists constructed by LGGa and LCGb respectively 
(the two initialising generators being given two different 
seeds). The TTT based estimate in table U is consistent 
with the table |TT] results of Deng and Blote most likely 
obtained from this generator. 

QTAa and QTAm are the QTA generator respectively 



TABLE I: Site percolation threshold estimates for the square 
lattice (pf(2048)) obtained by various pseudorandom number 
generators (PRNGs) as described in the text. 



PRNG 


Pf(2048) 


TTa 


0.59274627(11) 


TTHa 


0.59274588(11) 


TTTab 


0.59274628(12) 


SWBb 


0.59274617(17) 


QTAa 


0.59274588(17) 


QTAm 


0.59274603(17) 


QTBa 


0.59274610(17) 


QTBb 


0.59274621(17) 


XGx 


0.59274596(15) 


MTa 


0.59274593(17) 


MTm 


0.59274585(16) 


MTHm 


0.59274598(12) 


DMTmm 


0.59274597(08) 



initialised from LCGa and LCGm. Since the QTAa and 
QTAm results are consistent, there is no evidence that 
estimates from the QTA generator are especially sensitive 
to initialisation. The union of these two data sets gives 
an overall estimate of pf (2048) = 0.59274595(12) from 
the QTA generator. This value is consistent with earlier 
results, in table [TTl obtained by Ziff and Stell with this 
generator. 

QTBa and QTBb are the QTB generator initialised 
from LCGa and LCGb respectively. Since the QTBa and 
QTBb results are consistent, there is no evidence that es- 
timates from the QTB generator are especially sensitive 
to initialisation. The union of these two data sets gives an 
overall estimate of pf(2048) = 0.59274616(12) from the 
QTB generator. This is consistent with the TTa value, 
and also with Newman and Ziff's similar observation re- 
garding these two generators [2^, [2^ . Referring to table 
im the value is also consistent with the various threshold 
estimates obtained by Newman and Ziff using (at least 
in part) this generator. 

MTa is the MT generator initiahsed from LCGa. MTm 
is the MT generator initialised from LGGm. Since the 
MTa and MTm results are consistent, there is no evi- 
dence that estimates from the Mersenne twister genera- 
tor are sensitive to initialisation. The union of these two 
data sets gives an estimate of pf (2048) = 0.59274589(12). 
This is consistent with the MTHm result derived from 
the hashed generator in table [H and hence there is no 
evidence that correlations in the MT output sequence 
influence the measurement at this level of precision. 
Hence the MT generator appears to be an adequate 
choice for the current application. Further combining the 
MTHm data into the union gives an overall estimate of 
Pf(2048) = 0.59274593(8) from the MT generator. This 
MT result is inconsistent with those of the TTT and (un- 
hashed) TT generators. The difference in results with 
respect to the QTB generator is no more than could be 



TABLE II: Published estimates of the square site percolation threshold. The pseudorandom number generator(s) used are given 
where known. Generator T is a Tausworthe generator, while C is a congruential generator. TTT is the generator most likely 
used by Deng and Blote. References are provided for both the result and the generator whenever those come from different 
sources. Uncertainties are quoted as one standard deviation statistical errors, except in the semi-rigorous results of Balister, 
BoUobas and Walters (99.99% confidence bound) and of Riordan and Walters (99.9999% confidence bound). Only those results 
derived from currently accepted scaling relations are shown from the greater collection in Hu, Chen and Wu. This table is 
essentially a continuation of that appearing in Ziff and Sapoval [41|, there going back to 1960. 



Year 


Ref. 


Author(s) 


Method 


Generator (s) 


Resuh 


1986 


m 


Ziff and Sapoval 


Hull-gradient 


T 


0.592745(2) 


1988 


[11,32] 


Ziff and Stell 


Hull-gradient 


QTA 


0.5927460(5) 


1989 


[42] 


Yonezawa, Sakamoto and Hori 


Planar crossing 




0.5930(1) 


1992 


[33] 


Ziff 


Hull-crossing 


QTA 


0.5927460(5) 


1994 


[43] 


Hu 


Histogram Monte Carlo 




0.592(8) 


1995 


m 


Hu 


Histogram Monte Carlo 




0.5928(1) 


1996 


[45] 


Hu, Chen and Wu 


Histogram Monte Carlo 




0.59278(2) 


1996 


[45] 


Hu, Chen and Wu 


Histogram Monte Carlo 




0.59283(4) 


1996 


[45] 


Hu, Chen and Wu 


Histogram Monte Carlo 




0.59267(6) 


1996 


[45] 


Hu, Chen and Wu 


Histogram Monte Carlo 




0.5814(30) 


1996 


[45] 


Hu, Chen and Wu 


Histogram Monte Carlo 




0.6041(30) 


2000 


[28,29] 


Newman and Ziff 


Toroidal wrapping 


TT, QTB 


0.59274621(13) 


2000 


[28,29] 


Newman and Ziff 


Toroidal wrapping 


TT, QTB 


0.59274636(14) 


2000 


[28,29] 


Newman and Ziff 


Toroidal wrapping 


TT, QTB 


0.59274606(15) 


2000 


\28,29} 


Newman and Ziff 


Toroidal wrapping 


TT, QTB 


0.59274629(20) 


2000 


[28,40] 


Ziff 


Hull-gradient 


QTB 


0.5927465(2) 


2002 


m 


Ziff and Newman 


Planar crossing 


QTB 


0.5927464(5) 


2003 


[46,47] 


Martins and Plascak 


Toroidal wrapping 


C 


0.5927(1) 


2003 


[46,47] 


Martins and Plascak 


Toroidal wrapping 


C 


0.5929(3) 


2005 


[30,3^ 


Deng and Blote 


Cylindrical correlation 


TTT 


0.5927465(4) 


2005 


[30,3^ 


Deng and Blote 


Cylindrical correlation 


TTT 


0.5927466(6) 


2005 


[30,3^ 


Deng and Blote 


Cylindrical correlation 


TTT 


0.5927466(8) 


2005 


[30,3^ 


Deng and Blote 


Cylindrical correlation 


TTT 


0.5927468(10) 


2005 


[35] 


Balister, BoUobas and Walters 


Semi-rigorous 


MT 


0.5927(8) 


2007 


[36] 


Riordan and Walters 


Semi-rigorous 


MT 


0.59275(25) 


2007 


[37] 


Lee 


Planar crossing 


MT, DMT 


0.59274603(9) 



expected by chance in a data set of this size. The SWB, 
QTA and XG generator based estimates are consistent 
with that of the MT. 

DMTmm is the DMT generator with its two initial 
hsts independently constructed, each from one of a pair 
of seed words, by LCGm. The DMTmm result is con- 
sistent with the combined MT result, thereby indicating 
that any possible correlations between lower order bits in 
MT output words are insignificant at this level of preci- 
sion, or at least no worse than correlations in the higher 
order bits. This suggests that the single MT generator 
will be adequate for the purposes of this study. Com- 
bining all four Mersenne twister based data sets; MTa, 
MTm, MTHm, and DMTmm, produces an estimate of 
Pf (2048) = 0.59274595(6). This is consistent with the re- 
sult of Lee, in table [TTl obtained with this same mixture 
of generators but in the microcanonical approximation 
RL{n/N) K, i?L.„. The combined value does not alter 



any of the above conclusions regarding the consistency or 
otherwise of other generators with the Mersenne twister. 
Regarding the previous estimate of Lee, it was observed 
here that n/N , such that i?2048,n = 0.5 + k/L (inter- 
polating to non-integer n), usually exceeds pf(2048) by 
approximately 2 x 10~®. That being so, a revised esti- 
mate of the published result would be Pc = 0.5927460(1). 
This adjustment is much smaller than statistical uncer- 
tainties. 

Although the procedure used here differs from those of 
previous works, the results obtained are found to be con- 
sistent when the same pseudorandom number generators 
are used. However, given the use of a consistent method, 
it has been shown that the results thus obtained can dif- 
fer with the choice of generator. The level of PRNG 
sensitivity will be method dependent. The spread in re- 
sults seen here is not extreme as only reasonable quality 
generators have been used. 



The SWB, QTA, QTB_, 
backed by strong theory 



m 



XG and MT generators are 



and have been 
Ziff has per- 



17 




extensively tested elsewhere [17|, |2 
formed a sensitive hull generating walk test upon sev- 
eral generalised feedback shift register generators lrn|. 
Two-tap generators performed poorly in this test which 
concluded that they best be avoided for critical appli- 
cations. Certain quad-tap generators, particularly QTB, 
performed very well. Analysis indicated that QTB should 
outperform QTA in principle, although no obvious prob- 
lems were observed in the latter. The MT generator has 
passed tuned collision tests conducted by Tsang, Hui, 
Chow, Chong and Tso [24|. The LCGa generator failed 
those same tests. L'Ecuyer and Simard have recently 
performed thorough randomness tests upon a large as- 
sembly of PRNGs, including SWB, QTB, MT, XG and 
LCGa [Hi. The XG generator passed all tests, the MT 
failed in a very limited number of instances, QTB and 
SWB both failed a small number of times, and LCGa 
failed badly. TT was not specifically tested, although 
two-tap generators typically performed poorly. 

Results from the Mersenne twister generator have been 
consistent under different initialisation methods and ef- 
fective lattice enumerations (hash functions). With the 
observation that results from the Mersenne twister dif- 
fer from those of the two-tap lagged Fibonacci genera- 
tor, in the presence of evidence suggesting that the two- 
tap suffers from significant output correlations, and in 
the absence of evidence for any such correlations in the 
Mersenne twister output sequence, further Monte Carlo 
sampling within this exercise shall be performed exclu- 
sively with the MT19937 algorithm. Note that results 
from the SWB, QTA, QTB and XG generators are con- 
sistent with those of the MT. 



V. THRESHOLD DETERMINATION 

Having identified the Mersenne twister as a suitable 
PRNG for the problem, a more precise determination of 
the square site percolation threshold can now be made. 
This will be based upon Monte Carlo estimates of the 
microcanonical RL,n curves for 128 < L < 4096 (a span 
of some three orders of magnitude in N) . 

Data for L < 1024 was obtained exclusively from the 
MT generator initialised by LCGm, and Monte Carlo 
sampling was conducted with the algorithm of Lee [33] . 
Sampling domains were n e {8900 : 10500} on the 
L = 128 lattice, n G {37300 : 40400} on the L = 256 
lattice, n £ {152300 : 158500} on the L = 512 lattice, 
and n £ {615500 : 627600} on the L = 1024 lattice. 
The data at i = 2048 is the combined MTa, MTm, 
MTHm and DMTmm data from table [H As noted, that 
data was obtained with the same algorithm over the site 
occupation domain n G {2474000 : 2498300}. Due to 
hardware constraints, the L — 4096 data was obtained 
with the more memory efficient algorithm of Newman 
and Ziff [2^, [2^. For this algorithm, the entire domain, 



n G {0 : N}, is sampled, however observations were made 
only for n G {9920000 : 9969000}. Once again, the 
LCGm initialised MT generator was used. Lattices of 
L much more than 4096 could not be accommodated by 
the computer system used without substantial decreases 
in performance. Estimates at each L are based upon be- 
tween 1 X 10^ (at L = 4096) and 4 x 10^ (at L = 128) 
independent samples per occupation level, n. These re- 
quired the generation of between 10^^ (at L = 128) and 
10^^ (at L — 4096) pseudorandom numbers. 

As before, these microcanonical ensemble crossing 
probability curves, Rl,^ were transformed into canon- 
ical ensemble crossing probability functions, Rl{p), by 
the convolution of equation ([T]) . Because the various mi- 
crocanonical sampling domains all encompass ±12(TL.p 
of the convolution region about the critical point, the 
domain restriction induced error in Rl{p) is completely 
negligible for the values of p considered here. 

Several statistics were calculated from each Rl{p) 
curve. These were Ziff 's median-p critical point estimator 
|3a |. Pm{L), defined such that 



RUPn^iL)) = 1/2, 



(3) 



the Reynolds, Stanley and Klein real-space renormalisa- 
tion group cell-to-cell estimator [4^, |43|, Pcc{L,), defined 
such that 



Rl{Pcc{L)) ^ Rl/2{Pcc{L)), 



(4) 



the Ziff and Newman linear combination estimator [34| . 
Ph{L), defined as 



Ph(i) = {Pm{L) + apcc(i))/(l + a). 



(5) 



and the real-space renormalisation group cell-to-site fixed 
point estimator of Reynolds, Klein and Stanley [50], 
Pr{L), defined such that 



RL{p,iL))=p,{L). 



(6) 



Numerical estimates for these quantities are shown in 
table [ml 

The estimators Pm and Pcc are believed to approach 
their limiting values on the infinite lattice as L~^~^/'^, 
where v ~ 4/3 j33l . [34| . The estimator ph is believed to 
converge to its limit at a faster rate of 2^-i-"-i/'^^ where 
Ziff and Newman have determined a value oioj = 0.90(2) 
[34| for the scaling exponent proposed by Aharony and 
Hovi [5l|,[54]. The estimator p^ is believed to approach its 
limit as L~^''^ , a much slower rate of convergence than 
for the other estimators [3J]. Although each of these 
four limits is numerically equivalent to the percolation 
threshold pc, it will be useful to adopt a general notation 
indicating the origin of any threshold estimates. 

To second order, the finite size scaling relation for the 
median-p estimator is 



PniiL) «p* 



aL 



-l-l/iy 



bL- 



-Ijv 



(7) 



TABLE III: Site percolation threshold estimators on square lattices of various sizes L. Pm{L) is the median-p estimator, Pcc{L) 
is the cell-to-cell estimator, ph{L) is the linear combination estimator, and Pi{L) is the fixed point estimator. Results were 
obtained with the Mersenne twister pseudorandom number generator. 



L 



Pm{L) 



Pcc(i) 



PlAL) 



Pr{L) 



128 


0.59266108(21) 






0.59598352(23) 


256 


0.59272062(18) 


0.5928085(4) 


0.5927460(5) 


0.59467466(18) 


512 


0.59273860(15) 


0.5927651(4) 


0.5927462(4) 


0.59389258(15) 


1024 


0.59274377(14) 


0.5927514(4) 


0.5927460(4) 


0.59342699(15) 


2048 


0.59274528(06) 


0.5927475(2) 


0.5927459(2) 


0.59315051(06) 


4096 


0.59274573(10) 


0.5927464(3) 


0.5927459(3) 


0.59298626(10) 



FIG. 2: Parametrised fit of first order scaling theory (equation 
((7| with 6 = 0) to experimental data (pm(I') of table [lll| for 
the median-p critical point estimator. 



FIG. 3: Parametrised fit of first order scaling theory (equation 
([8]) with c = 0) to experimental data (pcc (L) of table IIII[) for 
the cell-to-cell critical point estimator. 
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A parametrised fit of equation ^ to the data of ta- 
ble Hn] produces p*^ = 0.59274595(4), a = 0.413(5), 
and b = 0.0(4). That the coefRcient b is indistinguish- 
able from zero suggests that a first order model (equa- 
tion ((T]) with 6 constrained to zero) is appropriate for 
the data. As such, a precise value for uj is unimpor- 
tant. In this first order case, the coefficients are evalu- 
ated as p*^ = 0.59274596(3), and a = 0.4135(7). The 
very good agreement between model and experiment is 
shown in figure [21 An empirical power law fit of the 
form pm(L) ~ p*^ - aL'' yields p^^ = 0.59274595(5), 
a = 0.42(2), and z = -1.75(8). That the value of z is in- 
distinguishable from the assumed exponent of —1 — 1/z^, 
further supports scaling of the form L~^~^^'^ as being the 
appropriate model for the data at this level of precision. 
All three p^ estimates are in good agreement with one 
another. 

The second order scaling relation for the ccU-to-cell 
estimator is given by 

pUL) « p:, + -L-i-i/"- + cL-i— V- (8) 

a 

where a = 1 - 2^1/'" [34]. The quality of the cell-to- 
cell data is lower than that of the median-p data, as 
each point is obtained from the intercept of two lines. 



with statistical uncertainties, at a shallow angle, and 
as Pcc{L) and PcciL/2) are not entirely independent. A 
parametrised fit of equation ^ to the data of table IIIII 
produces p*^ = 0.5927458(2), a = 0.441(5) and c = 
—9(8). Coefficient c is not inconsistent with zero, and a 
first order model (equation ([8]) with c constrained to zero) 
does fit the data, as shown in figure [3l with coefficients 
oip*^ = 0.5927459(2) and a = 0.417(4), in good agree- 
ment with the median-p estimator results. An empirical 
power law fit of the form Pcc(L) ^ p*^ — (a/a)L^ yields 
p*^ = 0.5927458(2), a = 0.34(8), and z = -1.71(4), con- 
sistent with the assumed L^^^^/" scaling relation. All 
three p*^. estimates are consistent with each other and 
with the estimates for pjj^, although the precision is sig- 
nificantly lower. 

The linear combination estimator of equation ([5|) was 
constructed by Ziff and Newman [3J| so as to cancel the 
first order terms of equations ([3]) and ([4]), leaving a faster 
approach to the percolation threshold; 



Ph(^) ^-^^1^ + 



1 



' T -I-Ul-I/V 



(9) 



(to first order). A parametrised fit of this expression to 
the data of table IIIII is shown in figure [5] and produces 



FIG. 4: Parametrised fit of scaling theory (equation ([9])) to 
experimental data {ph{L) of table [lll| for the linear combina- 
tion critical point estimator. 




FIG. 5: Parametrised fit of scaling theory (equation (|10|l ') to 
experimental data {pt{L) of table [III[I for the renormalisation 
group fixed point percolation threshold estimator. 
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pI = 0.59274596(7), and (b+ac) = 0.3(9). The threshold 
result is in good agreement with those obtained from the 
median-p estimator data. The value of (6 + ac) is also 
in agreement, although this is not saying much given the 
large uncertainties. That this value is essentially indis- 
tinguishable from zero is a reflection of the rapid rate of 
convergence of the Ph{L) estimator with L, as suggested 
by equation ^ , and the relative lack of precision in the 
PhiL) data. This is unsurprising given that no higher 
order terms were apparent in either the PmiL) or pcc{L) 
data sets. As such, the data was inadequate to empiri- 
cally test the assumed scaling exponent and is even con- 
sistent with ph{L) = constant = pjj, for which fitting the 
weighted mean gives pj^ — 0.5927460(1). 

The second order scaling relation for the fixed point 
renormalisation group estimator is given by 



p,{L)^p; + rL-^^-' + sL 



-2/1/ 



(10) 



[34| , and so has a slower rate of convergence to its infinite 
lattice limit than any of the other estimators above. A fit 
to the data of table IIIII yields well defined numerical val- 
ues for the coefficients; p; = 0.5927441(7), r = 0.1238(2), 
and s = —0.021(6). However, as shown in figure [5l the 
model of equation (jTU]) is but a loose match to the data at 
best, with higher order terms evidently remaining signif- 
icant. As such, the stated uncertainty in p* is misleading 
and will be addressed shortly within the next section. An 
empirical power law fit of the form p^{L) ^ p*—rL^ yields 
p* = 0.592743(2), r = 0.1219(8), and z = -0.748(2), 
consistent with the assumed first order exponent of —l/v. 
The first order model (equation ([9]) with s constrained to 
zero) returns p* = 0.5927456(8) and r = 0.1233(1). Of 
course, neither of these two functions describe the data 
any better than does the second order model. 



VI. ROBUSTNESS 

Several estimates have now been made for the square 
site percolation threshold, pc, using all the data of table 
IIIII and with varying degrees of precision. Of these, the 
most precise is pc — Pni — 0.59274596(3), obtained from 
the median-p estimator data using the first order scaling 
model pi„(i) = Pm~^0'L~^~^^'' ■ I* is prudent to establish 
the robustness of the results with respect to variations in 
the data and in the assumed model, over what domains 
the various models are valid, and how the domain and 
any fixed model parameters influence the estimate of Pc- 
There is a trade off between fitting to as much data over 
as great a domain as possible, so as to reduce statistical 
sampling fluctuations and hence to refine the result, and 
fitting to only data from large lattices where finite size 
effects are smaller and the scaling theories better describe 
the data. The results in tableUare consistent, where they 
overlap at L = 128 and L — 256, with those of Ziff and 
Newman [3J|. Hence their data was used to extend the 
domain down to L = 8 as necessary. 

The fixed point renormalisation group estimator is pos- 
sessed of good quality data, but has a slow rate of con- 
vergence to its limit p*. The Pv{L) data of table IIIII can 
be reasonably well fit with the addition of an L~^l^ term 
to the model, however the coefficient of L"''/'^, in an even 
higher order model, is not zero. Values of the coefficients 
fluctuate with the order of the model, suggesting that 
even higher order terms remain significant. Empirical 
power law fits are consistent with the leading order ex- 
ponent being — 1/i^, however a purely first order model 
does not fit the data well until the domain is truncated 
to L > 256. Results for p* are sensitive to the pres- 
ence or absence of individual data points, the L — 4096 
point altering the result by ±1 x 10~^. Under different 
models and data ranges, threshold estimates range from 
0.592744 to 0.592746. The difference is much larger than 
the uncertainty in the individual estimates and so not 



much weight should be given to those. Consequently, al- 
though the raw data at a given L is relatively precise, the 
slow rate of convergence of the fixed point renormalisa- 
tion group estimator leads to only a very rough figure of 
p* = 0.592745(f). 

The linear combination estimator suffers from rela- 
tively large statistical uncertainties in the data, and 
points are not entirely independent of one another. How- 
ever, the estimator does claim a very rapid rate of con- 
vergence to its limit, p^^. The model of equation ([9]) fits 
the data well for L > 32. Results thus obtained range 
from pI = 0.59274594(5) to p*^ = 0.59274603(8), with 
the presence or absence of individual data points mak- 
ing differences of as much as ±4 x fO~^ in pjj. Allow- 
ing for alternative values of the parameter w, between 
0.85 and 0.95, the estimate changes by no more than 
±f X f 0^*^. The data is not precise enough to either sup- 
port or falsify the assumed scaling relation, and is not 
inconsistent with ph{L) — constant. Even so, all esti- 
mates for pjj were consistent with one another and with 
the f28 < L < 4096 p\,{L) data mean of 0.5927460(f). 
ffence the linear combination estimator appears to be ro- 
bust, and the mean value, which covers the entire range 
of results, should be a more than safe estimate for pj^. 
The value oi pi = 0.59274596(7), obtained from aU the 
P\i{L) data of this study, should be reliable. 

With similarly low data quality, non-independent 
points, and a slower rate of convergence, the ccU-to- 
cell renormalisation group estimator should not be ex- 
pected to provide any refinement in pc over the lin- 
ear combination approach. Over domains where the 
various models fit the data, cell-to-cell results for p*^ 
range from 0.5927458(2) to 0.592746f (2). Sensitivity 
to the presence or absence of individual data points is 
as for the linear combination results, but here this is 
much smaller than statistical uncertainties. The estimate 
p*^ = 0.5927459(2), obtained earher from fitting the first 
order scaling model to all the Pcc data of table Ifffl is 
in agreement with the entire range of cell-to-cell results 
above, and so is robust, if imprecise. 

The median-p based estimates have the same rate of 
convergence as the cell-to-cell estimates, but with inde- 
pendent data points of much higher quality. The median- 
p estimates are less sensitive to the presence or absence 
of any one particular data point, this making a differ- 
ence of at most 2 x fO~®, and typically of less than 
f X fO~®, in the result for p^. The first order model 
fits the data for L > f 28, with results lying in the range 
p*, = 0.59274594(3) to p^ = 0.59274596(4). The second 
order model fits the data for L > f6, with results lying 
between p^ = 0.59274591(8) and p*^ = 0.59274600(5). 
The empirical power law model makes a good fit for 
L > 64, with estimates of p*i running from 0.59274589(2) 
up to 0.59274603(2), and scahng exponents in the range 

— 1.729(7) to —1.79(2). As noted in the previous section, 
the first order fit matches the data of table lllll verv well, 
the empirical fit agrees with the assumed exponent of 

— 1 — l/v, and coefficients of higher order terms were in- 



significant. This indicates that the first order model does 
indeed provide an accurate description for the finite-size 
scaling behaviour of the median-p estimator. The es- 
timate thus obtained, of p^ — 0.59274596(3), does not 
quite encompass the entire range of results above. Allow- 
ing for an extreme scenario, where even the model and 
scaling exponent may not be quite right, a more conser- 
vative figure of p^ = 0.59274596(4) does cover all of the 
above results. Hence this final value of the median-p es- 
timate for pc should be quite dependable. Incidentally, a 
standard error of 4 x f 0~^ in pc is approximately what 
would be expected from the total amount of data sam- 
pled in this study (as listed in the Pui{L) column of table 
ini[) . Parameter a of equation ([7]) shows much more sensi- 
tivity to the data domain and model than does p^. The 
fitted value given in the previous section was the most 
precise obtained. An overall result of a = 0.415(5) is 
more reasonable in light of the other estimates. 

The four estimators have now produced equally many 
robust estimates for the two-dimensional square site 
percolation threshold pc- As summarised in table IIVI 
these are p*, = 0.592745(1), p*^ = 0.5927459(2), pj; = 
0.59274596(7), and p*^ = 0.59274596(4), in good mutual 
agreement. Taking pc = 0.59274596, and returning to 
the canonical spanning probability curves, a good match 
between the data of 128 < L < 4096 and the theory of 
Rl{Pc) « 0.5 + kL-^ -t- 0(L-2) was had for k = 0.317(1). 
No higher order terms were seen, with the coefficient of 
L~^ being indistinguishable from zero. The value of k 
found here is a little lower than those of Ziff, k = 0.319(1) 
m, and Newman and Ziff, k = 0.320(1) [H. The dif- 
ference in pf(2048) resulting from using k = 0.3f7, as 
opposed to A; = 0.320, in equation ^ is around 6 x f 0^^. 
This is much less than the statistical uncertainties in the 
results of table [H upholding the claimed insensitivity of 
those estimates to k. Hence those results remain rea- 
sonable (PRNG biased) estimates for pc, and the direct 
comparison with earlier pc estimates is valid. The es- 
timate a = 0.415(5) found above is consistent with the 
results of Ziff, Newman, Hovi and Aharony [33, [SJ, [S^l 
(a here equates to their ratio bo/ai). Since k equates to 
bo, it follows that ai = 0.76(1) from the data obtained 
within this exercise. This estimate is also consistent with 
those of previous works ^, ^, ^] . 

The above results are based on data acquired solely 
from the Mersenne twister, that generator having been 
determined as suitable for this problem. In section 
HVl results obtained from the SWB, QTA, QTB and 
XG generators were found to be consistent with results 
obtained from the Mersenne twister. Although those 
four generators were not tested to the same extent as 
Mersenne twister, there is no objective reason to dis- 
count them entirely. Incorporating the data obtained 
from these generators earlier leads to revised values of 
Pm(2048) = 0.59274532(4), pr(2048) = 0.59315055(4), 
Pcc(2048) = 0.5927476(1), Pcc(4096) = 0.5927463(2), 
Ph(2048) = 0.5927459(1), and Ph(4096) = 0.5927459(3) 
for the various estimators of table IIIIl Note that the 
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FIG. 6: Parametrised fit of empirical power law model 
Pin{L) = p^ — aL^ to combined experimental data from the 
Mersenne twister, subtract with borrow, xorgens, and both 
quad-tap generators. 
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majority of the data remains Mersenne twister based. 

Use of these revised values does not alter either the 
fixed point limit, p*, or the cell-to-cell limit, p*^. The 
linear combination limit is raised to p^ = 0.59274598(6), 
an adjustment of rather less than its statistical uncer- 
tainty. 

A parametrised fit of equation ^ to the revised 
median-p data yields p*^ = 0.59274598(3), a = 0.415(7), 
and b — 0.1(5). As before, the coefficient of the higher 
order term is indistinguishable from zero. A first order 
fit (of equation ([7]) with b constrained to zero) produces 
p*^ = 0.59274598(3), and a = 0.414(1). An empiri- 
cal power law fit of the form Pm{L) ~ p^ — aL^ finds 
p*^ = 0.59274598(4), a = 0.41(2), and z = -1.75(1). 
The excellent agreement between this model and the ex- 
perimental data is shown in figure [6] The fitted value of 
z is indistinguishable from the assumed scaling exponent 
of —1 — l/i^ (with 1/ = 4/3). All fitted parameters are 
consistent across the three models. 

Performing robustness checks as before, the first or- 
der model fits the data for L > 128, with results ly- 
ing within a worst case range of p^ = 0.59274598(6) to 
p'^ = 0.59274599(8), and much more typically within 
p*, = 0.59274598(2) to p*^ = 0.59274599(4). The second 
order model fits the data for L > 16, with results lying 
between p^ = 0.59274596(3) and p^ = 0.59274602(4). 
The empirical power law model makes a good fit for L > 
128, with estimates oi p*^ running from 0.59274598(4) up 
to 0.59274599(8), and scaling exponents, z, in the range 
— 1.74(1) to —1.77(1). Hence the data supports the va- 
lidity of the first order model with the assumed scaling 
exponent, and a standard error of 3 x 10~® in p^ ap- 
pears fully justified. The various threshold estimates are 
summarised in table IIVI Using the revised data, and 
Pc — 0.59274598, the estimate of the finite size correc- 
tion parameter remains unchanged at k = 0.317(1). Nor 
is any significant change is seen in parameter a. 



TABLE IV: Infinite lattice limit estimates for the percolation 
threshold. Results are shown, by estimator, for the Mersenne 
twister only data (MT, MTH, DMT), and also for the com- 
bined generators data (MT, MTH, DMT, SWB, QTA, QTB, 
XG). 



Limit 


Mersenne 


Combined 


Pr 


0.592745(1) 


0.592745(1) 


Pec 


0.5927459(2) 


0.5927459(2) 


Ph 


0.59274596(7) 


0.59274598(6) 


Pm 


0.59274596(4) 


0.59274598(3) 



Assuming the suitability of the Mersenne twister 
PRNG for this particular Monte Carlo application, and 
also assuming that the median-p estimator approaches 
the critical point as Pm(i) — Pc esc L"^^^/"^, where 
v — 4/3, as supported by the data, then a robust es- 
timate for the square site percolation threshold is pc — 
0.59274596(4). A value for w is not required. Further as- 
suming the suitability of the subtract with borrow, xor- 
gens and both quad-tap generators, the additional data 
adjusts this estimate to Pc = 0.59274598(3). Continuing 
to assume the reliability of those generators, while drop- 
ping the assumed scaling exponent and requiring only 
that Pm(i) — Pc oc i^, for some z, the estimate becomes 
Pc = 0.59274598(4). These three estimates are mutually 
consistent to well within statistical uncertainties. The 
most precise of them has a standard error of 3 x 10~*, 
however a degree of caution is warranted in that none 
of the generators were tested to that level of precision. 
That being the case, this study's final estimate for the 
square site percolation threshold is 



Pc = 0.59274598(4). 



(11) 



This, primarily Mersenne twister based, estimate is 
consistent with almost all previous results in table |TI1 
In particular, it is in good agreement with the Mersenne 
twister derived estimate of Lee. Taken collectively how- 
ever, those results, excluding that of Lee, would suggest 
a higher value for pc, in the vicinity of 0.5927463(1). Al- 
though the value obtained here lies well outside of that 
range, the difference could be attributable to the vari- 
ous pseudorandom number generators used. While it is 
not impossible that the result obtained here may reflect 
some detectable influence of the chosen generators, pre- 
cautions against this were taken and no evidence of bias 
was found. 



VII. CONCLUSIONS 

Increasing availability of highly parallel computer facil- 
ities now makes it practical to obtain significant quanti- 
ties of Monte Carlo data from large lattices. This allows 
for greater precision in derived statistics, but requires 
very good quality pseudorandom number generators as 
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it is well established that inadequate generators lead to 
erroneous results. 

Tests were performed upon several generators and it 
was found that use of simple two-tap generators should 
probably be avoided for this application. The MT19937 
generator appeared to be suitable and was adopted for 
the majority of the Monte Carlo sampling conducted 
within this study. No dependence was found upon the 
(reasonable) choice of generator initialisation. 

Percolation threshold estimates subsequently made 
from various crossing probability statistics were found to 
be in good mutual agreement. The most precise of these 
was obtained from the median-p estimator. Data quality 
was such that precise results could be obtained without 
the need to assume a particular scaling exponent. Even 
so, results were in good agreement with a leading expo- 
nent of —1 — 1/v and no higher order term was found. 
The square site percolation threshold was subsequently 
determined to be Pc = 0.59274598(4). 

This estimate is consistent with the majority of earlier 
results on an individual basis, but not with those same 
results combined. Evidence suggests, however, that at 



least some of those earlier results have been influenced 
by the pseudorandom number generators used. The gen- 
erators used here appear to be of adequate quality, and 
the main generator, MT19937, passed an application spe- 
cific test of randomness. Furthermore, efforts were made 
to ensure the reliability of the error bounds in that final 
estimate, which should, then, be accurate. 
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